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Abstract 

Computation of the probability that a random graph is connected 
is a challenging problem, so it is natural to turn to approximations 
such as Monte Carlo methods. We describe sequential importance 
resampling and splitting algorithms for the estimation of these prob¬ 
abilities. The importance sampling steps of these algorithms involve 
identifying vertices that must be present in order for the random graph 
to be connected, and conditioning on the corresponding events. We 
provide numerical results demonstrating the effectiveness of the pro¬ 
posed algorithm. 


Notation 


G = {V,‘£) 

Graph with vertex set V and edge set “£. 

g(V) 

Subgraph of Q induced by the vertex set V. 

T(P) 

Set of subgraphs of Q induced by a vertex subset, or 
the power set of V. 

X 

The set of up vertices of a network. 

D, 

Vertices of g which are definitely up. 

Pr 

Vertices of g which are possibly up. 

x r ,x r 

Random subsets and supersets of X. 


Possible values of D r . 

K. 

Possible values of X r . 

F 

Event that X is connected. 

F r 

Event that X can still be connected conditional on D r . 


l 



I r (d r ) Indicator function of P (F | D,. = d r ) > 0. 

\Z\ The number of elements of set Z. 

B (x, r ) Closed ball of radius r around x. 

B (x, r ) Open ball of radius r around x. 


1 Introduction 


In its broadest sense, network reliability is the study of the performance 
characteristics of systems that can be modeled by random graphs. The 


most common application is the study of communication networks (Cancela 


et al 


Posse 

Lillo 



2009), but other applications include electricity networks (Chassin and 


2005 Pagani and Aiello 2013) and air transport networks (Zanin and 


2013; Wilkinson et al 2012; Cardillo et ah 2013). Often such systems 


are highly reliable, and the problem of estimating failure probabilities for 
these systems is one of estimating a rare-event probability. The most widely 
studied network reliability model is the K-terminal network reliability model , 
where the system is operational if a specified set of vertices is connected; see 


Gertsbakh and Shpungin (2010) for further details. 


We consider the problem of estimating the residual connectedness reliabil¬ 


ity (Lin and Ting, 2015), defined as follows. Let Q — (‘R, *E) be an undirected 
graph. Vertices of Q fail independently with probability 1 —p, and edges are 
considered failed if either of the vertices fails. The assumption that all ver¬ 
tices fail with equal probability is a notational convenience; generalization 
to the case where vertices fail with different probabilities is straightforward. 
Failed vertices are said to be down and functioning vertices are said to be 
up. The overall network is said to be UP if the up subgraph induced by 
the up vertices is connected, and DOWN otherwise. The probability that 
the system is in the UP state is called the residual connectedness reliability. 
Important applications include Radio Broadcast Networks and Mobile Ad hoc 
networks (Royer and Toh, |1999 


Tanenbaum 2002). 


In practice exact computation of the residual connectedness reliability 
(RCR) is found to be difficult. The formal statement is that computation of 
the RCR is NP-hard (Sutner et al 1991). That is, computing the RCR is 


at least as difficult as the hardest problem in the computational complexity 
class NP (Sipser, 1996). For NP-hard problems we generally turn to approx¬ 


imations such as Monte Carlo methods. 

Most established Monte Carlo methods for the related K-terminal network 
reliability problem , such as approximate zero-variance importance sampling 
fL’Ecuyer et al 2011), the merge process (Elperin et al 1991) and generalized 
splitting (Botev et al 2013) cannot easily be adapted to the RCR problem. 
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These methods generally rely on the network structure function being mono- 
tonic, which is not the case for RCR. One method that can be adapted is 


the Recursive Variance Reduction (RVR) method (Cancela and Urquhart 


2002). Another, specifically designed for the RCR model, is the conditional 


Monte Carlo method in Shah et al (2014). We propose a novel sequential im 


portance resampling algorithm for the estimation of the RCR, which can be 
interpreted as a “splitting” algorithm combined with importance sampling. 

The splitting method , first described in Kahn and Harris (1951), is a 
Monte Carlo technique for the estimation of probabilities of rare events. If 
the rare event can be written as the intersection of nested events, then the 
rare-event probability can be written as a product of conditional probabil¬ 
ities. These probabilities will typically be much larger than the rare-event 
probability and can therefore be estimated more accurately. The key fea¬ 
ture of the splitting method is that sample paths of some random process 
are replicated or split when they hit some “intermediate level of rareness”. 
Although the estimators of the conditional probabilities are generally not 
independent, under the splitting method their product is still an unbiased 
estimator of the rare-event probability. The splitting method and similar 
sequential importance resampling algorithms such as stochastic enumeration 


(Villcn-Altamirano and Villen-Altamirano 1991 i 

\u and Beck, 

2001 Liu 

et al 
2014 

2001; Del Moral et al, 2006 Rubinstein, 201 

0; Vaisman and Kroese, 

have found numerous applications including network reliability (Hector 

et al 

2008), queuing theory (Glasserman et al 1999 

Garvels, 200( 

)), particle 


Morio, 2013). 


transmission (Kahn and Harris, 1951) and air traffic control (Jacquemart and 


The rest of this paper is organized as follows. Section [2] introduces the 
RCR model and describes two existing Monte Carlo estimation methods for 
the RCR problem. Section [3] introduces two splitting algorithms for the 
estimation of the RCR. Section [4] gives an overview of the transfer matrix 
method, which allows the exact computation of the RCR in some cases. 
Section [5] describes a simulation study performed to compare the different 
estimation method and discusses the results. A list of proofs is given in 
Appendix [Aj 


2 The Residual Connectedness Reliability Model 

2.1 Definition 

Let Q = (T 1 , £) be a connected graph in which the state of a vertex v G 
V is a binary random variable. If the binary variable takes value 1 then 
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v is functioning correctly; otherwise the vertex has failed. For notational 
convenience we assume that these random variables are all independent and 
identically distributed. This means that each vertex is in the up state with 
probability p and the down state with probability 1 — p. The subset of V 
containing the up vertices is denoted by X. We will identify the set X with 
the subgraph Q (X) which has vertex set X and edge set 

{{D, v 2 } eT,-. v 1: v 2 E X} . 

Let 7 ('h') be the set of all subgraphs of Q induced by a vertex subset. Define 
ip as the binary function on 7 (‘V) which takes value 1 for connected graphs 
and 0 otherwise. When we write ip (X), this will be 1 if Q (X) is connected 
and 0 otherwise. We are interested in estimating the residual connectedness 
reliability 

l {Q,p) = P (X is connected) = P (ip (X) = 1). 

It will be convenient to take some enumeration V \,..., V\<p\ of the ver¬ 
tices V and use this as a total ordering. This makes it is possible to write 
statements such as V\ < v 2l and min {v \, v 2 ) = V\. We also use the nota¬ 
tion [v, oo) = {zu E 7* \ zv > vj. Note that the ordering used is completely 
arbitrary. 


2.2 Existing Methods 


The Recursive Variance Reduction (RVR) method of Cancela and Urquhart 


(2002) can be adapted to the RCR problem as follows. Let U be the event 

and let £ be the probability that the input 
The RVR method begins by writing 


that all vertices of Q are up 
random graph model is connected 


i = P (<p (X) = 1) = ip {Q) P (U) + P (<p (X) = 1 | U c ) P ( U c ) 

= p w +F((p(x) = i\u c ) (i-p^y (i) 


The event U c is the event that at least one vertex fails. Let V\ be the first 
vertex that fails (with respect to the total ordering of vertices). If no vertex 
fails then we say that V\ — oo, but U c occurring is equivalent to the condition 
Vi < oo. In this case we have 


P (Vi = Vi | Vi < oo) 
Equation (|TJ) can be rewritten as 


p l X (1 — p) 

1 — pID 


£ = p M + P (cp (X) = 1 | Vi < oo) (l - p 1 ^ 1 ) . 
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The problem now is to estimate P (tp (X) — 1 \ V\ < oo). But conditional 
on an observed value of V\, we know that vertex V\ is down , and all ver¬ 
tices smaller (with respect to the total ordering of vertices) than V\ are 
up. We have therefore fixed the states of some of the vertices. So condi¬ 
tional on Vi, we have a random graph where the number of vertices with un¬ 
known state is smaller than in the original random graph. Define £S l > (iq) = 
P(p(X) - 1 | Vy = ui). Then 


£ = + E [£ (1) {Vi) | Vl < oo] (l — . (2) 

The RVR method first simulates V\ \ V\ < oo. If £W (Vi) can be explicitly- 
computed, then an estimator for £ is 

p m + e m (vy(i-pM). 


If £^> (Vi) cannot be explicitly computed then the second down vertex V 2 is 
simulated and the expansion used in ([lj is applied recursively to £- 1 ’ (Vi). 
See Cancela and Urquhart (2002) for further details. 

Another approach is the simple conditional Monte Carlo method sug¬ 
gested in Shall et al (2014). In this method a random order is selected for 
the vertices. The states of the vertices are then simulated in that order, until 
the first up vertex is identified, denoted by c o. The connected component 
for the identified up vertex is then simulated. At this point the states of 
some vertices are still unknown. The probability that the random graph is 
connected is the probability that the connected component of oj is the only 
connected component. This is the probability that all vertices with unknown 
state are in fact down, which is a power of 1 — p. See Shah et al (2014) for 
further details. We will refer to this method as the conditional Monte Carlo 
method. 


3 Splitting for Residual Connectedness Reli¬ 
ability 

We construct a fixed-length Markov chain, denoted by {D,.}f =0 . Increasing 
values of r correspond to increasing information about the random graph X. 
At the last time stage we have complete information about X, so in fact X 
is equal in distribution to D^. In some cases we will be able to determine 
that (p (X) = 0 before “time” R. We will use this Markov chain to construct 
a classical splitting algorithm. Note the use of r for the time variable; this 
is intended to emphasize that our “time” variable corresponds to the radius 
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(a) Representation of X as a sub¬ 
graph of §. Solid dots represent up 
vertices. 

Figure 1: Example value of X (left), 
graph Q (right). 
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(b) Ordering of vertices of 

Ordering of vertices for the 5x5 grid 


of a certain disk. As in Section [2] we will assume some arbitrary ordering of 
the vertex set ©’. 


3.1 Notation 


We illustrate the ideas in this section using the 5x5 grid graph, with R = 3 
and the value of X shown in Figure l(a)l Recall that the value of X will be 


the last value of the Markov Chain (DJto- We assume the vertex ordering 


shown in Figure 1(b) 


We first construct a metric on the set ©' of vertices. For vertices Vj and 
Vj, the distance d(Vi,Vj) is defined as one less than the number of vertices in 
the shortest path in Q from v % to v r For example, the shortest path from i) % 
to Vi is {Vi}, so d (Vi , Vj) = 0. If Vj and v 3 are connected by an edge then the 
shortest path is {v i: Vj}, so that d (v t , v 3 ) — 1. As we assumed that Q was 
connected, such a path always exists. 

Let R be some positive integer. We now define random sets X”, for 
integers 0 < r < R and n > 1. Recall that X is the set of up vertices, and 
define X), = min X. That is, X) is the set containing the smallest up vertex 
with respect to the partial orde ring (which does not depend on r). For the 
value of X given in Figure 1(a), X,, 1 = {2}. For a subset x of vertices, define 


Next 


= min {x G X | d i 


> R — r,x > max x} . 


In words, Next (x, r) is the next up vertex, after those in x, which is more 
than distance R — r away from the vertices in x. Such a vertex may not exist, 
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(f) Value of X 3 . 


Figure 2: Vertex subsets with R = 3 and the value of X shown in Figure 
|l(a)[ The red crosshatched vertices are up vertices contained in the specified 
vertex subset. Solid vertices represent other vertices which will eventually 
be determined to be up. 

in which case Next (x, r) = 0. We now make the inductive definition 

X” +1 = X" u Next (X?, r). 

This means that X” is an increasing sequence of sets of up vertices, each set 
containing at most n vertices. 


Example values of Xq,Xq an d Xj) are shown in Figures 2(a) - 2(c) As 


Q is a finite graph, at some point this sequence must reach some largest 
set, and stop increasing. We therefore define X r = U^LiX”- Note that if 
r = R then in fact X# = X. In words, X r is the set generated by taking the 
first up vertex, and then iteratively selecting the next (w.r.t. ordering) up 
vertex that is at least distance R — r away from those previously selected, 
until no such vertex can be found. Example values of X 0 , X], X 2 and X 3 


are shown in Figures 2(c) - 2(f). The set of possible values of X r will be 


denoted by 3£ r . We view the generation of X r as the result of a function 
GenerateSubset (X, R, r ) defined in Algorithm [I] 
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Algorithm 1: GenerateSubset (X, R, r) 
input : Set XCf, maximum radius R, radius r 

output: X r — ^ 

1 if X = 0 then 

2 return 0 

3 X, 4 0 

4 while 3s G X such that s > maxX r and d (s, X r ) > R — r do 

5 |_ X r 4— X r U min {s G X | s > maxX r , d (s, X r ) > R — r} 

6 return X r 


The notation is intended to emphasize that the random sets {X r } are 
random subsets of X. The {X r } are all strongly dependent on X and on each 
other. The intermediate random sets {X"} are only needed to construct X r , 
and we rarely refer to them past this point. 

For every possible value x r of X r let 

Up r (x r ) = |J (B(x,R-r) n [x, oo)) . (3) 

i£x r 

The set Up r (x r ) represents all vertices of Q that could possibly be up given 
that X r = x r . Define 


X, = Up r (X r ). (4) 

Note that the random sets {X r } are random supersets of X (see Proposition 


A.l). The value of X r completely determines the value of X r because of the 
functional relationship in (J4j) . As X r and X r are subsets and supersets of X 


we have that 


X r C X C X r , 


for all 0 < r < R. 


( 5 ) 


Finally, we take partial intersections and unions to define the supersets 
and subsets 


D r = (J X*) and P r = f|X t . (6) 

t =o t=o 

The mnemonics D and P are chosen because conditional on known Xi,..., X r , 
the set D r represents the set of vertices that are definitely up, and P r repre¬ 
sents the set of vertices that are possibly up. The set of possible values of D, r 
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(b) 

Value of D 
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Figure 3: Values of Di and D 2 with R = 3 and the value of X shown in 
Figure 1(a) The red crosshatched vertices are up vertices contained in the 
specified vertex subset. Solid vertices represent other vertices which will 
eventually be determined to be up. 


will be denoted by 3> r . Example values of Di and D 2 are shown in Figure 
|3(a)| and 3(b) Note that the value of D 0 is always equal to Xo which is 
shown in Figure 2(c). The value of D/j is always equal to X, which is shown 


in Figure l(a^ 


ft is clear from (|5]) that X r = x r implies that x, r C X C Up r (x, r ). 
However from Proposition |A.5| the converse is also true. So 

P (X, = x r ) = P (x r CXCUp r (x r )) = P (x r C X, Up r (x r ) c C X c ) 

= pix r i n _ p yu Pr (x r ) c i _ 

From ([5]) and (J6]) we see that D r and P r are subsets and supersets, so that 
D r CXC P r , 0 < r < R. (7) 

For every 0 < s < r we know that 

X s C D, C X C P r C Up, (XJ . (8) 


Proposition A.5 and (J8J) imply that knowledge of D r completely determines 
the value of X, for 0 < s < r. Specifically, X s = GenerateSubset (D r , i?, s ). 
Consequently D r also completely determines the values of D s and P s for 
0 < s < r. 

This makes (D,}f=o a Markov chain, as its state at any time completely 
determines its sample path up until that time. By contrast, {X r }f =0 is not 


9 






















a Markov chain. For example, consider Figure 2(c) and 2(d) Information 
from X 0 shows that vertex 14 of Q is up, but this is not reflected in X,. So 
if X 14 is the binary variable controlling the hnal state of vertex 14, we have 


P (X 14 = 0 | X 0 = x 0 ) = P ( X 14 = 0 | x 0 = x 0 , Xi = Hi) = 1, 
P(X 14 = 0|X 1 =x 1 )^l. 


This means the Markov property fails to hold. 


3.2 Splitting Algorithm 

Let F be the event that X is connected. Let F r be the event that X can still 
be connected, conditional on the observed value of D r . More precisely, 

F = MX) = 1}, 

F r = (3x G 7 (F) with D r C x C P r , such that ip (x) = 1} 

= (P (F | D r ) > 0} . 

Note that the subsets D, are increasing and the supersets P r are decreasing, 
so that {-F r }f =0 i s a decreasing sequence of events, with F R = F. 

For every 0 < r < R and d r G define I r (d r ) = I {P (F | D r = d r ) > 0}. 
That is, I r is the indicator function that takes value 1 if X can still be con¬ 
nected given the observed value of D, , and 0 otherwise. Note that we can 
evaluate I r (D r ) by performing a depth-first search of P r , started from only 
one of the vertices of D r . Let x DF s be the vertex set traversed by the depth 
first search. Then xdfs is a connected graph, and if xdfs contains D r we 
have D, C xdfs Q P, This makes xdfs a connected graph which has a 
non-zero probability of occurring, in which case 


P (F | D r = d r ) > P (X = xdfs I Dr = d r ) > 0. 
The nested events {F r }^ ={d allow us to write P (F) as 

P(F) =P(Fi | F 0 )P(F 2 | Fi) • ■ • P (F r | F r _0. 


We can therefore estimate P (F) as a product of estimates of conditional 
probabilities using a splitting algorithm. This leads to the fixed splitting al¬ 
gorithm given in Algorithm [2] For notational simplicity we take the splitting 
factors to be integers. Removing this restriction is conceptually easy, see 
Section 10.6 of Kroese et al (2011) for further details. 
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Algorithm 2: Splitting algorithm for RCR 
input : Connected graph Cj. maximum radius R. reliability p. integer 
splitting factors k 0 ,... , hn_\. initial number of particles N 
output: Estimate of RCR 

1 * 3 / i — 0 // Samples that have hit an intermediate level 

2 for i = 1 to N do 

3 sample D from Do 

4 if Iq (D) then add D to 


5 

6 

7 

8 
9 

10 


for r = 1 to R do 

3? i — 0 // Samples that have hit the next level 

for d r _i G W do 

for i = 1 to fc r -l do // Split each particle k r -\ times 

sample D from (D r | D r _i = d r _i) 
if I r (D) then add D to J 


n 


swap and 2? 


12 return \ < 3/\ ^A r n^=o' 


3.3 Importance Sampling and Resampling Algorithms 


Experimentally we observe that most particles reach event Fr- 1 , but then 
tend not to reach event Fr = F. That is, we have decomposed the rare- 
event probability into a product of R conditional probabilities, where the 
last conditional probability is very small and the others are relatively large. 
Typically in splitting we would solve this by adding events between Fr_i and 
Fr. But due to the discrete nature of the constructed Markov chain it is not 
obvious how this can be done. Fortunately, the particles which reach Fr -i 
have a very special structure. 

In graph theory a cut vertex or articulation vertex is a vertex of a graph 
which, when removed, increases the number of connected components. Now, 
consider the set V cut of cut vertices of Xr_i. If Fr -i occurs and some v G 
Kut is n °t contained in X R _,, then v must be up if X# = X is to be 
connected. For a proof of this statement see Proposition |A.6| Note that 
this only applies at step R — 1, and not at any other step. Figure [| gives a 
graphical representation of this statement. Figure 4(a) shows P 2 , which is 
the set all of vertices that are possibly in X. Figure 4(b) shows the smaller 
set D 2 , which is the set of all vertices that are definitely in X. Figure |4(cJ 
shows all the cut vertices of P 2 , and Figure 4(d)] shows the single cut vertex 
of P 2 which is not in D 2 . If this vertex is down then X cannot be connected. 
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(c) Cut vertices of P 2 . 
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(d) Cut vertices of P 2 that 
are not in D 2 . 


Figure 4: Visual representation of the result described in Proposition A.6 


with R = 3. Figure [4(d) shows the single cut vertex of P 2 that is not in D 2 . 
This vertex must be up in order for X to be connected. 
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(a) The single biconnected 
component of the value of P 2 
shown in Figure |4(a)| which 
contains more than two ver¬ 
tices. 


(b) Another biconnected 
component of the value of 
P 2 shown in Figure [4 (a) | 


Figure 5: Biconnected components of the value of P 2 shown in Figure 4(a) 
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Algorithm 3: Sequential importance sampling algorithm for RCR 
input : Connected graph Cj. maximum radius R. reliability p. integer 
splitting factors k 0 ,..., hn_\. initial number of particles N 
output: Estimate of RCR 

1 *3/ i — 0 // Samples that have hit an intermediate level 

2 for i = 1 to N do 

3 sample D from Do 

4 if / 0 (D) then add D to 


5 

6 

7 

8 
9 

10 


for r = 1 to R — 1 do 

2f i — 0 // Samples that have hit the next level 

for d r _i G W do 

for i = 1 to k r —i do // Split each particle k r -\ times 

sample D from (D r | D r _i = d r _i) 
if I r (D) then add D to J 


n 

12 

13 

14 

15 

16 

17 

18 


swap and 2? 

"W i — 0 // Importance weights 

for d fi _! G & do 

cut cut vertices of Up ij ._ 1 (y) 

for i = 1 to kft-i do // Resample each particle ku-i times 

sample D ~ (D^ | Dr_i = d^-ijCut C X) // imp. sampling 
if Ir(D) then add pl cut \ d «- i l to W] // Add imp. weight 

return (rir^Lo 1 ) ^iwGW w ^ Return average of imp. weights 


Conditioning on these cut vertices all being up leads to the importance 
resampling algorithm shown given in Algorithm [3j Note that the condition¬ 
ing step is a type of importance sampling and leads to an algorithm involving 
weighted particles. The factor of pR cut \yl in the final estimator is the impor¬ 
tance weight. 

The cut vertices of a graph can be used to decompose it into biconnected 
components. A biconnected graph is a graph which has no cut vertices, and 
therefore cannot be disconnected by the removal of a single vertex. A bi¬ 
connected component is a maximal biconnected subgraph of some underlying 
graph. The value of P 2 shown in Figure 4(a) has only one biconnected com¬ 


ponent with more than two vertices. This subgraph is shown in Figure 5(a^ 


and another component is shown in Figure 5(b) Note that the biconnected 


components are not disjoint and share vertices. Any vertex contained in more 
than one biconnected component must be a cut vertex. 
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Proposition A. 7 shows that if Fr -i occurs, then X is connected if and 
only if V cnt C X and the ‘restriction’ of X to every biconnected component 

of P R _i is a connected graph. But the connectivity of X restricted to one 

component is independent of the connectivity of X restricted to another 


component. We can use this independence to estimate P (F 
more efficiently. This leads to Algorithm [4] (SIS). 


Dfl-i — d/j_i 


Algorithm 4: Sequential imp. sampling algorithm for RCR (SIS) 
input : Connected graph Cj. maximum radius R, reliability p, integer 
splitting factors k 0 ,..., k R _i, initial number of particles N 
output: Estimate of RCR 

1 Identical to lines 1 - 14 of Algorithm [ 3 ] 

2 for d.R_i G ^ do 


3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


cut f— cut vertices of Up ij ._ 1 (d#_i) 

bi •?— biconnected components of Up R _ 1 (d^_i) 

count! 0,..., count|bi| f— 0, resid f— 0 

for i = 1 to |_ b (/do // Simulate most copies efficiently 
sample D from (Dr | Dr_i = d^_i,cut C X) 

for j — 1 to |bi| do 

| if (p (D fl bij) then county <— county + 1 

for i = [ |b ^_iJ |bl1 + 1 to k R _i do // Simulate remaining copies 
sample D from (Dr | Dr_i = d^_i, cut C X) 

if ip (D) then resid <— resid + 1 

|cut\d H _ 1 | f j-j-M county _|_ resid^ to W 


add p' ( 


14 return [ N 


-1 


E 




W 


We can also apply the same conditioning ideas at steps other than R— 1. 
This has the additional complication that not every cut vertex of X r is in 
fact required to be up in order for X to be connected. However if t; cut 
splits X r into components, at least two of which contain vertices of X r , then 
u cut must be up. This property can be checked by performing a depth-first 
search started at all vertices adjacent to u cut . This much more ‘aggressive’ 
conditioning leads to sample paths with significantly different weights. It 
is therefore sensible to substitute resampling for splitting, as this tends to 
automatically remove sample paths with small weights. These ideas lead to 
Algorithm [ 5 ] (SIR). 
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Algorithm 5: Sequential imp. resampling algorithm for RCR (SIR) 
input : Connected graph Q, maximum radius R , reliability p, initial 
number of particles N 
output: Estimate of RCR 

1 i — 0 // Samples that have hit an intermediate level 

2 'W i — 0 // Importance weights 

3 for i = 1 to N do 

4 sample D from Do 

5 if I 0 (D) then add D to add 1 to W 


6 

7 

8 
9 

10 

11 

12 

13 

14 

15 


for r = 1 to R do 

averageWeight «- \W\~ X E„, e r 

Sf <— with repl. sample of size N from (&, 'W") // Resampling 

W <- 0 , «- 0 

for d r _! £ f do 

cut G- cut vertices of Up r (d r _!), cond 0 
for c G cut do 

components connected components of Up r (d r _i) \ c 
n |{comp G components: comp D d r _i 7 ^ 0}| 
if n > 2 then add c to cond 


16 

17 

18 


sample D from (D,. | D r _! = z,cond C X) // imp. sampling 
if I r (D) then 

| add D to add averageWeight pl cond \ z l to W 


19 return N 1 YLw&r w 


// Return average of importance weights 


4 The Transfer Matrix Method 


In order to compare the different Monte Carlo methods in Section [5] it is 
useful to use graphs that are reasonably large and for which the RCR is 
exactly known. Assume that for 0 < i < |'R| the number of vertex subsets 
with 1 vertices that induce a connected subgraph is q. Then 


M 


PMx) = i) = £p(qx) = i||x|=; 

i =0 

Rl 




v l (1 - p) 


M-i 


X 

i =0 


c i m 


1^1 


p l (i-py 


(1 - p ) 


M-i 


i =0 


Computation of the values q has been shown to be complete (Sutner 


et al 1991). However for certain classes of regular graphs the transfer matrix 
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method from statistical physics (Klein et al 1986 Kloczkowski and Jernigan 


1998; Clisby and Jensen 2012) can be used to compute the integers c* rela¬ 
tively quickly. Once the numbers c % are known, the RCR can be computed 
quickly for any value of p, to arbitrary accuracy. 

As an illustration of the transfer matrix method assume that we are 
interested only in computing the total number of connected subgraphs of 
the 5x5 grid graph. The graph can be considered as five vertical ‘slices’ 
consisting of five vertices each. As we move from left to right the ‘state’ of 
the graph at a particular point consists of the states of those five vertices 
and information about the paths that connect them which lie to the left. Six 
of the 52 possible ‘states’ are shown in Figure |6j In the first example state 
there are three vertices present, which are all connected together by a path 
lying to the left. In the fourth example state there are three vertices present, 
and the first and last are connected by a path lying to the left. The middle 
vertex is not connected to the other two by a path lying to the left. In the 
third example state all three vertices are connected by a path which lies to 
the left, although the lower two vertices are necessarily connected as they are 
adjacent. These 52 states include two ‘empty’ states where no vertices are 
present; one which occurs before any ‘state’ which contain a vertex, and one 
which occurs after a ‘state’ which contains a vertex. 

When a connected subgraph is decomposed in this manner there are only 
a certain number of possibilities for the initial and final states. For example 
the first state in Figure [6] is not possible as an initial state and state 4 as a 
final state would result in a disconnected subgraph. Further only a limited 
number of transitions between states are allowed. The number of connected 
subgraphs can therefore be written as x T B 4 y. Here x is a binary vector 
containing the value 1 if a state is a possible initial state for a connected sub¬ 
graph, and y is a binary vector containing the value 1 if a state is a possible 
final state for a connected graph. The 52 x 52 binary matrix B contains the 
value 1 if a transition between two states is possible for a connected subgraph. 
This approach can be applied in general to regular planar graphs. By keeping 
a running count of the number of vertices we can determine the number of 
connected subgraphs for every total number of vertices. Our implementation 
made use of the Eigen linear algebra library (Guennebaud, Jacob et al 2010), 
the MPFR numeric library (iFousse et al, 2007) and Boost C++ libraries. 


5 Numerical Experiments 


We performed a simulation study to compare four Monte Carlo (MC) meth¬ 


ods for the RCR problem. The four methods were conditional MC (Shah 
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Figure 6: A sample of the six of the 52 possible states for the transfer matrix, 
for the 5x5 grid graph. 


et al 2014), the RVR method of Cancela and Urquhart| ( [2002| , the sequential 
importance sampling algorithm given as Algorithm|4|(SIS) and the sequential 
importance resampling algorithm given as Algorithm [5] (SIR). 

We previously assumed that the reliability of every vertex was equal. Not 
only is this notationally convenient, but it allows us to compute the exact 
RCR using the transfer matrix method if Q is a grid graph. The largest 
graph for which we applied the transfer matrix method was the 11 x 11 
grid graph. Although there are 2 m possible subgraphs, the transfer matrix 
method allows the computation of the c* in around an hour. The largest 
number of connected subgraphs occurred when i = 75, for which the number 
of connected subgraphs was on the order of 10 31 . Note that this is well 
outside the maximum integer value of 2 64 — 1 representable on a computer 
without more specialized software. 

We compared the four methods on square grid graphs with length 8, 9,10 
and 11. The parameter p was varied between 0.1 and 0.65 in steps of 0.05. 
At larger and smaller values of p the target event was no longer rare and 
we do not consider these cases interesting. In addition we also applied all 
methods to the 14 x 14 grid graph, for which it is not possible to compute 
the unreliability exactly. 

The methods were compared on the basis of their estimated relative error 
(RE) and work normalized relative variance (WNRV). For an estimator p of 
a probability p which takes on average time T to compute, these quantities 
are defined as 


RE (/)) = 


A/Var (p) 


and 


WNRV (p) = 7 


p p* 

In order to estimate the RE and WNRV, each method was run 100 times. 
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Figure 7: Relative error results for the 8x8 grid graph. Shading represents 
bootstrapped 95% confidence intervals. 


Sample sizes for every method were 10 6 per run. For SIS the splitting factors 
are a required input to the algorithm. A separate run with sample size 
10 6 and splitting factors identically 1 was performed in order to estimate the 
splitting factors. These estimated values were then used for all the subsequent 
100 runs. The time required to estimate these splitting factors is not included 
in the work normalized results. 

The simulation study shows that the RVR method is not competitive with 
the other three methods tested. For the 8x8 grid graph the average result 
from the RVR method is numerically accurate for the exactly computed value, 
although the estimated RE (Figure [T]) and WNRV (Figure [8]) show that the 
RVR method is inferior to the other methods. The shading in Figure [7] shows 
95% confidence intervals obtained by bias-corrected and accelerated (BCA) 
bootstrapping (Davison and Hinkley |l997) using the R package boot (Canty 


and Ripley, 2014). In the case of the 11 x 11 grid graph the RVR method 


also performs poorly, with the added problem that for p between 0.2 and 
0.35 the RVR method is not close to the exact probability (Figure [9]). 

In every case either SIR or conditional MC has the lowest relative error. 
Which of these performs best depends on p. In general for p < 0.25 con¬ 
ditional MC performs best but for p > 0.25 SIR performs best. A similar 
pattern is observed for the work normalized relative variance, although SIS 
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Figure 8: Work normalized relative variance 8x8 grid graph. 



Figure 9: The absolute empirical bias as a proportion of the true value, over 
100 replications for the 11 x 11 grid graph. 
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and SIR and conditional MC are somewhat closer in terms of performance. 
The value of p — 0.25 separating these behaviors is believed to be an artifact 
of onr choice of grid values for p. The value of 0.25 is the value that is closest 
to the value p* of p that solves 


p — E 


|X|M 


-1 



That is, the threshold value p* is the value of p for which the expected 
proportion of up vertices under the conditional distribution is exactly p. For 
example, Figure 10 shows the expected proportion of up vertices for the 
11 x 11 grid graph, with a vertical line at p* = 0.2454. Note that in Figure[9j 
for p < p* SIS and SIR do not always converge numerically to the true value 
using 100 replications, even on a log scale. Similarly, for p > p* conditional 
MC does not always converge to the true value using 100 replications. 



Figure 10: Expected proportion of up vertices for the 11 x 11 grid graph, 
under the conditional distribution. The blue line represents y = x, the red 
line represents p = 0.2454. 

The simulation results suggest that very different strategies must be em¬ 
ployed to estimate the RCR, depending on the value of p. At some value of 
p* (which we assume is unique) the expected fraction of up vertices under the 
conditional distribution is equal to p*. For values of p slightly smaller than 
p* the expected fraction of up vertices is much smaller than p, and for values 
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of p slightly larger than p* the expected fraction of up vertices is much larger 
than p. This occurs because for p smaller than p* the ‘easiest’ way to ob¬ 
tain a connected graph is to delete all but one of the connected components 
obtained under the unconditional measure. However if p is large then the 
easiest way to obtain a connected graph is to add extra vertices, connecting 
up the components generated under the unconditional measure into a single 
component. SIS and SIR take the approach of adding extra vertices, and are 
therefore only useful when p > p*. SIR is more aggressive than SIS in adding 
vertices, and the difference in performance between the two cases is therefore 
more pronounced. Conditional MC can be interpreted as ‘deleting’ vertices, 
and is therefore efficient for p < p*. 

Figure [TT] gives the relative errors for the four tested algorithms for the 
11 x 11 grid graph. The shaded regions indicate 95% confidence intervals 
for the RE for the SIR and conditional MC methods, obtained using BCA 
bootstrapping. The relative errors of SIR and conditional MC differ by a 
factor of up to 100, and which method performs better depends on whether 
p < p* or p > p*. On a work normalized basis (Figure [l2j) conditional MC 
outperforms SIR by a factor of up to 10 5 for p < p*, but SIR outperforms 
conditional MC by a factor of up to 10 2 for p > p*. 



H 


0.1 ' ok ' 0.3 ' 0.1 ' 0.5 ' 0.6 


P 


Figure 11: Relative error results for the 11x11 grid graph. Shading represents 
bootstrapped 95% confidence intervals for Conditional and SIR. 
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Figure 12: Work normalized relative variance for the 11 x 11 grid graph. 


Figure [13] shows the relative errors of the tested algorithms for the 14 x 14 
grid graph. The shaded regions indicate 95% confidence intervals for the RE 
for the SIR and conditional MC algorithms, obtained using BCA bootstrap¬ 
ping. We do not show results from the RVR method as the estimated relative 
errors were orders of magnitude higher than those for the conditional, SIR 
and SIS algorithms. For the purposes of calculating the relative errors we 
used the estimates given by conditional MC for p < 0.25 and the estimates 
given by the SIR method for p > 0.25. 

Note that for p < 0.25 the SIR method appears to have lower relative 
error. However the relative error estimates are believed to be inaccurate in 
this case, similar to the situation for the RVR method on the 11 x 11 grid 
graph (Figure [TTj) . The true RE for the SIR method is believed to be orders 
of magnitude larger than the RE for conditional MC for p < 0.25. The 
opposite situation occurs for p > 0.25; the relative error of conditional MC 
appears to decrease as p increases to 0.35, but this is believed to be because 
the RE is poorly estimated in this case. The true RE for conditional MC is 
believed to be orders of magnitude higher than the RE for the SIR method, 
for 0.3 < p < 0.45. The relative errors are believed to be well estimated for 
all three methods for 0.45 < p < 0.6, and in this region the relative error of 
the SIR method is a factor of 100 smaller than the RE for conditional MC. 
For the cases where the relative errors are poorly estimated, it is believed to 
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be computationally infeasible to run the relevant algorithm enough times to 
obtain an accurate estimate. 



Figure 13: Relative error results for the 14 x 14 grid graph. Shading repre¬ 
sents bootstrapped 95% confidence intervals for the SIR and conditional MC 
algorithms. 


Note that the SIR and SIS algorithms are not well suited to the situation 
where Q is a large subset of a complete graph. In that case every vertex is 
adjacent to every other vertex. So D 0 ,.... D#_i represent knowledge of only 
the first up vertex, and represents knowledge of the state of every vertex. 
In this case the SIR and SIS algorithms we describe are similar to crude MC 
or the RVR method. Further work is required to find methods that behave 
well in this situation. Another direction for further work is to find algorithms 
that efficiently estimate the RCR in the regime where p is small. Although 
conditional MC works acceptably well, it is somewhat unsophisticated. For 
the special value p = p* (which appears to be around 0.25 in our examples) 
none of the algorithms we tested work well. As shown in Figure 11, the SIS, 
SIR and conditional MC algorithms all had a relative error close to 1 in that 
case. This case appears to be particularly difficult as methods based around 
increasing or decreasing the number of up vertices will not be effective. 
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A Proofs 

Proposition A.l. The random variable X r is a superset ofX. 

Proof. Recall that X r Q X r . Assume that there is a vertex x' of X which 
is not in X r and therefore not in X r . If %' > maxX r , then we must have 
d(x',X r ) > R — r to avoid having x' G X r . So Next (X r , r) = x'. This 
implies that x' G X r , which is a contradiction. 

Assume that x' < maxX r , and define the set 

N = {n G Z + | x' > maxX"} . 

There is a maximum element l of this set, implying that x' > maxX ? r and 
x' < maxX^ +1 . As x' is not in X r we must have d(x',X r ) > R — r. This 
means that Next (X^ r ) = x' , in which case we would have x' G X r , which 
is a contradiction. 

Therefore every vertex of X must be contained in X r . □ 

Proposition A. 2. Take any 0 < r < R and x r G Stf r . Then for any x so 
that x r C x C Up r (x r ), it holds that minx r = minx. 

Proof. Clearly minx r = minUp r (x r ) by the definition of Up r in (|3]). But 
x, r C x C Up r (x r ) implies that 

minx r > minx > minUp r (x r ) . 


So minx r = minx. □ 

Proposition A.3. Assume that x = {xi,... ,x n } is an ordered sequence of 
vertices so that d (xk, {xi,..., Xfc_i}) > R — r. Then f%f_ r consists exactly of 
all such sequences x. 

Proof. It is clear that GenerateSubset (x, R, r, =, x) for any such x. But the 
definition of GenerateSubset in Algorithm [l] will always generate vertices 
with the stated property. □ 

Corollary A.4. 7/x G S£ T , then the value obtained by removing any number 
of points is still in 3P r . 
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Proposition A.5. Take any 0 < r < R and x r G Then for any x so 

that x, r C x C Up r (x r ), 


GenerateSubset (x, R, r) = x r . 

Proof. Let y = GenerateSubset (x, R , r) and let ay,..., x Ux be an enumer¬ 
ation of x r according to the arbitrary ordering. Similarly let y 1 , ... ,y ny be 
an enumeration of the set y. From Proposition A.2 we know that y 1 — x r . 
All the subsequent points in y are more than distance R — r from y 1 — x{. 
Let A = (y 1 , oo) = (x 1 , oo). Then 


y fl A C x C Up r (y D A) , X r PliCxC Up r (X r fl A). 


Taking the intersection with A removes only the first point of both y and X r , 
so by Corollary A.4 these are still in SF r . Applying Proposition |A.2 again 
shows that y 1 = x 2 . Applying this argument inductively shows that n y = n x 
and y = x . □ 


Proposition A.6. Assume that d#_i is a value of Dr_i which implies that 
F r _i occurs and let p R -i be the associated value for P R _i. If v is a cut 
vertex o/pr_i then 


F Fl — d^_!} C — d/j.x} fl {v G X} . 


Proof. If v G D/j_i the statement is trivial as D/?_i C X. So assume that 
v g D^_i. 

As v is a cut vertex, it separates pr_i into two non-empty components, 
denoted by C\ and C 2 - Assume that component C\ does not contain any 
vertices of d^_! and let zv be a vertex in Cf. By the definition of p_r_i we 
must have zv G xr_i. This implies that there is a vertex of x R _x contained 
in B ( zv , 1) C C\ U v. As we assumed that v fL x^^, we must have a vertex 
of Xr_! contained in C\. 

This implies that C\ contains a vertex of d^_x, and by an identical argu¬ 
ment so does Ci- So there are vertices of X in both C\ and C 2 , and if v fL X 
these vertices will lie in different components of X. □ 

Proposition A.7. Assume that d R _! is a value of Dr_! for which F R _i 
occurs and let p R -i be the associated value for Pr_i. Let V CVLt be the set 
of cut vertices of Pr-x, and let B = be the set of all biconnected 

components of p R -i- Define the events 

|B| 

A = {Dr.! = dR_!} , B = {Kut C X} , C = f| {ip (X n Bf) = 1} . 

i =1 

The notation Xfl Bi refers to the subgraph of X induced by the vertex set B j. 
Then F nA = AnB nC. 
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Proof. Let H be the graph with vertices B, and edges between Bi and Bj if 
B t fi Bj ^ 0. Note that as Pr_i was connected, H must be too. 

Assume now that events A, B and C all occur. Take any vertices B j 
and Bj of H which are connected by an edge, and pick arbitrary vertices 
Vi G Bi fl X and Vj G Bj fl X. Let c be some vertex contained in Bi D Bj. 
Note that the event B implies that in fact Bi fl Bj C X, as all vertices in the 
intersection are cut vertices. 

As Bi fl X is a connected graph there exists a path P t in Bi fl X from 
Vi to c and another path Pj in Bj fl X from c to Vj. This gives a path in 
X connecting v, and Vj. The connectivity of H means this argument can be 
extended to the case where Bi and Bj are not adjacent in H. This proves 
that, given our assumptions, X is connected. We have therefore proved that 
AnBnCcFnA. 

Now assume that F and A occur. Proposition |A.6| implies that F fl A C 
Ap\B, so event B must occur. Now pick any biconnected component Bi and 
any two vertices V\ , v-2 G Bi fl X. As X is connected there exists a simple 
path between V\ and v- 2 . If there is such a simple path P that leaves Bi fl X 
then the subgraph PU Bi of P^-i is a biconnected component strictly larger 
than Bi. This would contradict Bi being maximal, so P must stay within 
Bi fl X, making Bi fl X connected. Therefore event C must also occur. 

This proves that F P\ A = An B P\ C. □ 
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